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We find that the hypothesis made by Jan, Stauffer and Moseley [Theory in Biosc, 119, 166 (2000)] 
for the evolution of sex, namely a strategy devised to escape extinction due to too many deleterious 
mutations, is sufficient but not necessary for the successful evolution of a steady state population of 
sexual individuals within a finite population. Simply allowing for a finite probability for conversion 
to sex in each generation also gives rise to a stable sexual population, in the presence of an upper 
limit on the number of deleterious mutations per individual. For large values of this probability, 
we find a phase transition to an intermittent, multi-stable regime. On the other hand, in the limit 
of extremely slow drive, another transition takes place to a different steady state distribution, with 
fewer deleterious mutations within the asexual population. 
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I. INTRODUCTION 

In a previous paper Q we have shown that a finite 
steady state sexual population may arise from a purely 
asexual one, if an excess of deleterious mutations causes 
the individual to engage in sex as a means of escaping 
death ^]. Since sexual reproduction led, in our model, 
to diploidy, it implied also the adoption of a mechanism 
for preferential expression of certain genes, and we as- 
sumed deleterious mutations to be recessive. Under 
various different assumptions regarding the subsequent 
mode of reproduction (i.e., whether sexual reproduction 
is hereditary or not) and of the number of offspring, we 
found that the diploid population always persisted, and 
that it was consistently more successful in escaping the 
effects of deleterious mutations. 

In the case where sexual reproduction was only prac- 
ticed as a means of escaping death from too many dele- 
terious mutations, but diploid cells were also allowed to 
multiply by mitosis, diploid individuals completely took 
over the population. Thus, in one of our models (Model 
I) |lj , we were able to demonstrate a possible scenario for 
the evolution of the analogue of a 'haploid - diploid cycle" 
Q where the organisms produce by asexual reproduc- 
tion as long as they are reasonably fit (or the conditions 
reasonably favorable) but engage in sexual reproduction 
when the going gets tough. || [| 

In this paper we will test whether a threshold mecha- 
nism for switching to sexual reproduction is necessary for 
the successful establishment of a sexual population. We 
simulate two strategies for the evolution of sex within a 
fixed population N of simple organisms, who are all ini- 
tially asexual (and haploid), and subject to a constant 
rate T of random mutations. Both haploid and diploid 
organisms die when the number of their expressed dele- 
terious mutations exceed a certain number. (Henceforth 
we will omit "deleterious" where it is evident that we 
mean a change away from the wild type gene.) 



The first strategy (Model ^4) is the adoption, with a 
certain probability, of sexual reproduction and conse- 
quent diploidy when the number of mutations exceeds 
a threshold, threatening extinction. The second strat- 
egy (Model B) involves a small but constant probability 
a for the accidental conversion to sex, independently of 
the number of mutations (or, equivalently, the fitness) 
of the individual. The cloning of sexual individuals is 
not allowed in cither Model A or B. On the other hand, 
we have tested for the effect of hereditary (habitual) v.s. 
non-hereditary (non-habitual) sex. 

In the implementation of both models, we have 
adopted a more realistic set of rules for the mechanism 
of dominance, that is, the expression of mutated alleles, 
than in our previous paper jlj . Here we allow a mutated 
gene to be expressed if the cell is homozygous for mutated 
alleles at this locus. Hence, the number of expressed dele- 
terious mutations for diploid individuals is the number m 
of different loci at which the cell is homozygous for mu- 
tated alleles. 

We find that both strategies A and B lead to a finite 
steady state sexual population, with typically a smaller 
average number of mutations (greater fitness) than the 
asexual population. Thus no threshold mechanism seems 
to be necessary for a successful sexual population to take 
hold. However, for habitual practice of sexual reproduc- 
tion by diploid individuals (i.e., those that are not facing 
extinction in Model ^4) calls for unrealistically large mu- 
tation rates in order for a macroscopic sexual population 
to be established in the steady state. 

The organisation of the paper is as follows. In the 
next section we explain in detail the two models and we 
report the result of our simulations. In Section III, wc 
display and examine the mean field evolution equations 
and discuss our findings in the light of these equations. 
In Section IV, we investigate the limits of strong and 
extremely weak driving of this system, for T — > 1 and 
r — > 0, as well as a transition to chaos via an intermit- 
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tent route, found for large values of a. A discussion of 
the results from similar models and directions for further 
research are provided in section V. 



II. MODELS FOR CONVERSION TO SEX 

We represent the genetic code of each one-celled indi- 
vidual with a bit-string of "0"s and "l"s, after the Penna 
model 0. At each locus, we have taken the value "0" 
to correspond to the wild type and "1," to a deleterious 
mutation (which we will call "mutation," for short, where 
this is not liable to lead to any confusion.) We use the bit 
defining the "sign", to specify whether the individual is 
asexual (+) or sexual (-). For asexual, haploid, cells, we 
have one 15-bit string, whereas, for the sexual cells, we 
have two 15-bit strings which are allowed to be different, 
i.e., the individuals are now diploids. 

A mutation consists of flipping a randomly chosen bit 
except the sign bit, and it is implemented by scanning all 
the individuals in the population, and, with probability 
r picking those to be mutated. Clearly there may be 
any number of mutated individuals at any one generation 
(time step) , the number fluctuating around T N, where TV 
denotes the total population. 

The number of deleterious mutations to is simply the 
number of "l"s for a haploid individual. For a diploid, 
the number of "expressed" deleterious mutations is taken 
to be the number of loci at which both homologous alle- 
les are set to "1." This is how we model the mechanism 
of dominance of the wild type (or, equivalently, the re- 
cessiveness of deleterious mutations.) We will use the 
term "fitness" loosely, for L — m; thus increasing to will 
decrease the fitness of the individual. 

In the steady state, the distribution of the asexual and 
sexual populations over to, are independent of T, for T > 
1/N. The cases where T < and T « 1 have interesting 
consequences, and are discussed in section IV. 

The probability of survival as a function of to is given 
by a Fermi-like distribution 0, P(m), 



P(m) 



1 



exp[/3(™ — /i)] + 1 



(1) 



For large /3 (or "low temperatures," in the language of 
statistical mechanics), P(m) behaves like a step func- 
tion [0. Individuals with to > /i die, those with m < fi 
survive, and those with m — /i survive with a probability 
of 1/2. In the simulations we set (3 = 10 and [i — 4. 

We keep the total population constant, as in the Red- 
field model H , by making up for the deficit in the popu- 
lation after all the bacteria have been either found fit for 
survival or killed off according to the survival probability 
in Eq. (|l]). Asexual individuals multiply by simply mak- 
ing another copy of themselves, namely by mitosis, while 
a pair of sexual organisms each contribute one bit-string 
to their offspring and die in the process. 

We performed the simulations on a fixed population of 
N = 1000, for 16-bit strings. The total number of time 



steps in each simulation is taken to be much larger than 
the time necessary for the transients to die off and the 
system to settle down to a steady state. Since the prob- 
ability to mutate a single gene in a diploid individual is 
r/(2i), on the average the steady state is reached after 
2L/T time steps, where 2L is the total number of genes 
in a diploid individual, or, in other words, the number 
of mutated genes in the population is greater than the 
total number of genes of one individual. In all the sim- 
ulations, the reported results are averages over 10 runs. 
The fluctuations depend on the model chosen, however 
the relative error estimate based on one standard devia- 
tion is typically less than about 6%, as long as there is 
only one fixed point for the dynamics. We will report on 
situations where we encounter an intermittent route to 
chaos in Section IV. 



A. Asexual steady state 

We start with a set of N initially identical asex- 
ual individuals, all identical to the wildtype, i.e., all 
0's. Under the conditions outlined above, without in- 
troducing sex, the population of asexuals settles down 
to the steady state distribution given in Q| for T > 
namely, n H [m)/N A = 0.012,0.098,0.356,0.531,0.001 for 
to = 0, . . . , 4. In this region this steady state distribu- 
tion is independent of T, which only sets the scale of time. 
That this should be so, is not self-evident, and only fol- 
lows from the form of the solution to the set of evolution 
equations, as shown in Section III. 



B. Triggering sex 

The alteration of the sex gene can be accomplished in 
two different ways. One can choose to trigger sex with 
a threshold mechanism or define a constant probability 
for each individual to become sexual. These mechanisms 
are further discussed in the following subsections. In ei- 
ther case, the haploid organism first makes a copy of its 
own set of genes, as if it were going to perform mito- 
sis, but then forms two gametes instead. One of these 
gametes will pair up with a gamete from another indi- 
vidual who has been turned on to sex, and the other will 
be discarded. One should note that sexual reproduction 
may be implemented in different ways, resulting in dif- 
ferent numbers of offsprings produced Q. In this paper, 
we will define sexual reproduction in such a way that 
when two sexual individuals mate they always give rise 
to one sexual offspring; thus, the population is reduced 
by one, each time an act of sexual reproduction takes 
place. Clearly, increasing the number of offspring will in- 
crease the advantage that the sexual population enjoys. 
Indeed, judging from our previous results pj], the number 
of offspring exceeding two would lead to the takeover of 
the population by the dipoid sexual types. 

When two diploid cells engage in sexual reproduction, 
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they each contribute one gamete towards a single diploid 
sexual offspring. Let us denote the two gametes as {Aa} 
in one parent, and {Bb} in the other parent. Then the 
genome of the offspring may be, {AB}, {Ab}, {aB} or 
{ab}. We do not allow for crossover between the gametes 
during sexual reproduction. 



We see that for T ~ 10 /N the proportion of the sexuals 
in the population saturates to ~ 70% as shown in Figure 
1, and remains at this value independently of the value 
of r. In order to obtain points near r ~ one has to 
do very long runs to get accurate results, and these are 
discussed in Section IV, as well as the chaotic behaviour 
displayed when T becomes too close to 1. 



1. Sex at the threshold of extinction - Model A 

In model A, alteration of the sex gene takes place only 
under special conditions, namely the threat of death due 
to too many mutations g| . Once the asexual steady state 
is reached, we allow the sex gene to be "turned on" for 
the least fit members of the population. In any pass 
through the population, if those individuals that are in 
the tail of the distribution (i.e. those with m > \x mu- 
tations) survive, then they are turned sexual by deter- 
ministically and irreversibly switching their sign bits to 
one. Once their sex bit is turned on, these individuals 
will be labelled "sexually active" and mate with other 
sexually active individuals. If there is only one active 
sexual at a certain time step then it must wait subse- 
quent generations until it finds a partner. After mating, 
the sexual individual becomes sexually inactive and the 
only way for it to become sexually active again is to face 
extinction once more. The deficit in the population due 
to deaths and to sexual reproduction is then made up by 
copying randomly selected asexual individuals. 
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FIG. 2: The distribution of both sexuals and asexuals over 
the number of expressed deleterious mutations m, for Model 
A. F = 10~ 3 . Hereditary sexuality is not allowed and the 
distributions are normalized to unity over each population 
separately. 
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FIG. 1: The percentage of sexual population v.s. T is plotted 
for Model A where hereditary sex is not allowed, for a popu- 
lation of 1000 individuals. The inset shows a larger range of 
r where the step-function like jump is more apparent. Both 
curves represent averages over 10 runs. 

In this model, therefore, there is no hereditary sexual- 
ity: there is, however, a hereditary transition to diploidy. 
This gives an unfair advantage to the sexuals in that they 
both enjoy the benefits of diploidy and escape the disad- 
vantage of 2 — > 1 reproduction. 



FIG. 3: The percentage of the sexual population v.s. F for 
Model A with hereditary sexuality introduced. The total pop- 
ulation is 1000 individuals and the results are averaged over 
10 runs. 

The steady state distributions of both asexuals and 
sexuals with respect to m are also independent of T, (See 
Figure 2) for T > j? and sufficiently smaller than 1. The 
peak of the distribution shifts towards lower m values for 
sexuals as a result of the salutary effect of dominance 
in diploidy, and the reshuffling effect of sexual reproduc- 
tion§. 
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Model A with Hereditary Sex 

We have also tested the case of hereditary, or habitual, 
sex, in which sexually active individuals can mate ran- 
domly either with sexually active individuals who have 
been converted to sex in that generation, or with individ- 
uals who have already been converted in some previous 
generation. As in the case of non-hereditary sex, the 
population is allowed to grow back to its fixed value by 
cloning randomly selected asexual units. 

This small difference results in a noticeable increase in 
the number of matings at each time step, and therefore 
leads to a decrease in the number of sexual individu- 
als in the steady state. We have found that the steady 
state comprises a macroscopic sexual population only for 
T > 1/N. For T < 1/N, the average number of sexual 
individuals drops to about 1%, or around 10 individu- 
als in a population of N = 1000. The sexual popula- 
tion increases linearly with T and reaches only ~ 15% 
(as compared to 70% for non-hereditary sex) as T ~ 1 
(see Fig. 3). The m-distributions are shown in Figs. 
(4a, b) for the asexual and sexual populations. The peak 
of the sexual population has shifted to 1 as a result of 
the greater number of mating events. Thus we may con- 
clude that hereditary and habitual sex in this model is 
punished more severely; the relative improvement in the 
mean value of m does not compensate sufficiently for the 
loss of the parents. 
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FIG. 4: The distribution of the a) asexual and b) sexual pop- 
ulation with respect to m, for different values of T for Model 
A with hereditary sexuality. The histograms are normalized 
to unity. 



with other sexually active individuals of that generation. 
(If there is only one active sexual at a certain time step 
then it has to wait till it finds a partner at a subsequent 
generation.) If we take sexual reproduction to be non- 
hereditary, after mating the sexual individual becomes 
sexually inactive. (Within some subsequent generation 
it can once more become sexually active with probability 
cr). The deficit in the population is made up by copying 
randomly selected asexual individuals. 
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2. Mutating the sex gene with constant probability - Model 
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Our second strategy for conversion to sex involves a 
constant probability a for the accidental conversion to 
sex, independently of the distance, as expressed by m, 
from the wildtype. For this model (Model B), once the 
asexual steady state is reached, at each generation we 
allow the sex gene to be "turned on" irreversibly, with 
a small probability a for each individual. Like in Model 
A, these individuals will be "sexually active" and mate 



FIG. 5: Percentage of the sexual population v.s. a/T plots 
for various F values for Model B. Hereditary sexuality is not 
allowed. All the points collapse onto a single curve in the 
interval shown. 

We find, (see Fig. 5) that this scenario again gives 
rise to a steady state macroscopic population of sexuals 
- but it is smaller than the one in Model A. The total 
percentage of sexuals is a function of a/T, as can be seen 
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from the figure, and grows with a/T. In Fig. 6(a,b), we 
display the distribution of asexual and sexual individuals 
over the effective number of mutations to, for two small 
values of a and T. The characteristic sandpile like Jl(J 
distribution of asexuals is accompanied by a distribution 
of sexuals which is again shifted towards smaller values of 
to. It is interesting to observe that raising a increases the 
total number of sexuals, and therefore depresses the num- 
ber of asexuals, as is to be expected. However, it is not 
immediately obvious why keeping a fixed and decreasing 
the overall mutation rate should decrease the number of 
asexuals. Clearly, raising T increases the death rate of 
both types of organisms, but since the conversion to sex 
is not coupled to the increase in the number of mutations, 
an increased T only benefits the asexuals who get cloned 
to make up the deficit population. For large values of a, 
a novel phase transition takes place, which is the subject 
of Section IV. 
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FIG. 6: Distribution with respect to m, for both sexual and 
asexual populations for Model B, without hereditary sex. a) 
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Model B with Hereditary Sex 

If the conversion to sexual reproduction is hereditary, 
then at any given time step all the sexual individuals 
mate, except for the odd guy out. In Fig. 7 we show the 
total percentage of the sexual population as a function 
of a alone. One sees that the growth is very close to 
linear with a, however the collapse as a function of a/T 
does not occur here. The curves extrapolate to zero at 
(7 = 0. As long as a > 1/N one may have a small but 
nonvanishing sexual population. For smaller values of a, 
the number of sexual individuals again fluctuates very 
strongly and is of 0(1). (see Section IV). 



III. MEAN FIELD EVOLUTION EQUATIONS 

To try to understand analytically some of the features 
found from the simulations, we have examined the be- 
haviour of the iterative equations that can be obtained 
for the different densities involved. 

Let us define the mutation matrix for haploids, 
T m , m+1 (T) = r(i-m)/iandT mim _ 1 (r) = Tm/L. Note 
that X^5=±i ^m,m+i(r) == r. All the other elements of 
this matrix are zero. 

For low temperatures and for fi, the upper limit of the 
number of mutations tolerated by the haploid individual, 
being set to four, the survival probability is given by, 



P(to) 




, TO = 

0, m > 4 



(2) 



A. Asexual steady state 

The time-evolution equations for the asexual popula- 
tion, with njy(TO-) being the number of individuals with 




FIG. 7: The percentage of the sexual population v.s. a for 
various F values for Model B with hereditary sex. The growth 
with a is linear for the different F values. 
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m mutated genes, are 

n H (m, t + 1) = (1 — r)rifl-(m) 

+ ^ T m+s . m n H (m + 5,t)-[l- P{m)]n H (m, t) 
a=±i 

+ JZt 1 -•PK)]»a("* / ,*)nfl(m,t)/JV > i . (3) 

The last term is the source term, arising from the replace- 
ment of the deceased individuals by randomly cloning the 
extant ones and Na = ^2 m n H (m) is the total number of 
asexual individuals. 

For large (3, one effectively has, 

n H (m,t+l) = (l-T)n H (m,t) 

+ ^2 T m+s , m n H {m + S,t) 
S=±l 

+ -n H {A,t)n H {m,t)/N A , (4) 

for m < 4. The source term [1 — P(4)]n# (4)nu(m)/NA 
has been replaced by its value \n^ (4)n_f/(m)/A r A, and 
it is assumed that nnirn > 4) = 0. This assumption is 
supported by numerical data in the steady state. 

Note that for FN ~ 0(1)), n H {A) will be small, i.e., 
of the order of unity. For m — 4, this enables us to put 
the source term in the last equation equal to zero, since 
it will be of 0(1/N) while the other terms are of O(l), 
and we get, 

n H (4,t+l) = (l-T)n H (4,t) 

+ T4+8.4n H (4: + S,t)- Trnff (4, t) 

4=±l 1 
= (l-r)n ff (4,t)+r(l-3/L)n ff (3) 

- \n H {^t) . (5) 

Then we see that in the steady state, one may re- 
place n#(4)/2 appearing in the source terms by — 
3/X)nir(3) — n.jj(4)]. This leads to equations that are ho- 
mogenous in r in the steady state, yielding, therefore, a 
steady state distribution of the population between sex- 
ual v.s. asexual individuals which are independent of T 
at least for T > 1/N. (see Fig. 1) Iterating these equa- 
tions leads to a steady state with an m-distribution that 
is in agreement with the simulation results 0. 



B. Coexisting asexual and sexual populations 

We now define a new quantity, Ud(to) as the number 
of m-mutation strings that belong to a diploid organism. 
The expected number of diploid organisms with m ex- 
pressed deleterious mutations can be obtained, once the 
71d(to) are known. 

The probability for two strings with mi and m 2 muta- 
tions (i.e., bit set to "1") to give rise to m loci at which 
both bits are "1" can easily be calculated. It is given by 



p{m; mi, ma) = 

mi!?«2!(£ — mi)!(i — m.2)! 
L\m\(m\ — m)!(m2 — m)\{L — mi — TO2 + m)! ' 

(6) 

for L—mi —m2 + m > and otherwise. This expression 
is symmterical in mi and m%, both of which must be 
> m. The number of diploid organisms with m expressed 
mutations is then, 

^ L L* 

n s {m) = - Y Y P (m;m 1 ,m 2 ) 

x n D (m 1 )n D (m 2 )/(2N s ) , 

(7) 

where Ns is the number of diploid organisms, 
Sm=o n s(™), and L* — min[L, L + m — mi]. The factor 
of h out front comes from converting from the number 
of gametes that are members of diploid organisms with 
m expressed mutations, to the number of such diploid 
organisms. The factor nD{rn2)/(2Ns) in the sum is the 
probability of encountering a gamete with m% mutations 
as the other member of the pair making up the diploid 
organism. 

A similar computation leads to the number of diploid 
individuals who die as a result of too many mutations, 

^ L L L* 

D D = ^Y1 Y X! t 1 ~ P(m)]p(m;m 1 ,m 2 ) 

rn— mi-m 777-2 —rn 

x n D (m 1 )n D (m 2 )/(2N s ) , (8) 

where L* us defined as above. 

The number of gametes with m mutations, which get 
removed because they happen to be members of diploid 
organisms which die, is 

L min[m,m"] 

dm = ^ ^ [1 - P(m')]p(m';m, m") 

m" — m'—O 

x n D (m,t)n D (m",t)/(2N s ) • (9) 

We must also define the number of gametes with m 
bits set to "1," that can take part in sexual reproduction, 
which is 

L 

dm = ^2p(i;m,m / )n D (m,t)n D (m',t)/(2N s ) (10) 

where L = min[L, L + 4 — m\. Since d m is only defined 
for m > 4, L = L + 4-m. Note that J2,n=4 = 2n s (4). 

Here we have only considered the scenarios without 
habitual sex. 

Model A 

We now have, from (H,||), for sufficiently large /3 
njj(m, t + 1) = (1 — T)n H (m, t) 
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+ ^ T m+ s, m n H {m + S,t)- S m ^nn(i, t) 
s=±i 

+ [|t»h(4, t) + D D (t) + in a (4, i)]»ff(n». . 

(11) 

The terms proportional to T are due to random mu- 
tation. The coefficient of the Kroenecker delta 5 m ^ is 
Uh (4) since all of the asexuals with m = 4 are removed 
either due to death or conversion to sexuals. The fi- 
nal term represents the number of m-mutation haploids 
which get cloned to keep the population constant; the 
expression in the square brackets is the number of in- 
dividuals which get removed from the population and 
determines the strength of this source term. The (3/4) 
factor multiplying n#(4) comes from two parts: (1/2) of 
the haploids with 4 mutations die; the other half is con- 
verted to sex, and mate, their number being once more 
halved as a result, contributing (l/4)n#(4) to the "re- 
movals." Djj (which is = (l/2)n a (4) for large 0) is the 
number of diploids that die, and (l/4)n s (4) comes from 
half of the m = 4 diploid population being converted 
to sex, their number being once more halved when they 
mate. 

The dynamics of the number of strands no(m) that 
make up diploid organisms is, 

1. 



n D (m, t + 1) = (1 - -T)n D (m, t) 

S=±l 

- d m (t) - -d, 



S,t) 



5 mA P(A)n H (A,t) (12) 



For diploids, the probability of a mutation hitting any 
one gene is halved, because there are twice as many of 
them. The d m term is the number of m-gametes that 
are removed as a result of death, and in practice (for 
large (3) is nonzero only for m > 4. The next term gives 
the reduction in the number of m-gametes as a result 
of sexual reproduction. A factor of (1/2) comes from 
the probability to engage in sex, and another from the 
fraction of gametes that are discarded as a result. Finally, 
there is a contribution from the conversion of haploids to 
diploids. We have neglected the situation where a) there 
is only one active sexual individual is present, so that 
no mating with concomitant discarding of a gamete, can 
take place; or b)& conversion from haploid to diploid is 
impeded because there is only one haploid strand with 4 
mutations. It can be checked explicitly that Eqs.([Tl 12) 
conserve the total population. 

Iterating these equations leads to a steady state dis- 
tribution that is roughly comperable but not identical to 
the simulation results (see Table I). For T = 10~ 3 the 
percentage of the sexual population is 24% of the total, 
and saturates to 36% as T is increased, as compared to 
70% from the simulations. This discrepancy seems to 
come from the fact that the dynamics is really driven by 
the strongly fluctuating small population at m = 4, and 
mean field theory is simply not able to capture this. 



TABLE I: The distribution of the population with respect to 
the number of expressed mutations, obtained from an itera- 
tion of the mean field equations for Model A. 
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0.0 
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Total 
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The distribution over m is also modified; one sees that 
the distribution of the asexuals is quite similar to the sim- 
ulation results, while the peak of the sexual distribution 
has shifted to m = 1, from m = 2. This indicates that 
the mean field theory overestimates the effect of remix- 
ing, as is to be expected, since the gametes, instead of 
being paired in a definite way at any given moment, are 
perpetually part of a single gene pool. 

Model B 

In this case we have a uniform probability for conver- 
sion to sex. The equations become, 

n#(m, t + 1) = (1 — r)rifl-(m) 

+ ^ T m+s , m n H (m + 6) - [1 - P(m)]n H (m, t) 
s=±i 

- an H (m) + {^[1 - P{m')]n H {m') + -<xN A 



D D (t) + -aNs(t)}n H (m)/N A 



(13) 



Here, haploids are converted to diploids and removed 
at the rate of cr, and the reduction in the population 
due to mating of recent converts gives the \(tNa term 
in the source. The sexuals moreover mate among each 
other with probability er, which leads to a further sink 
with strength ^aNg. Apart from these, the terms are 
identical to Eq. (JTT|) . The dynamics of the m-gametes 
are, 

n D (m,t+ 1) = (1 - ir)n D (m,t) 



^ T m+s , m {-T)n H (m + S,t) - d r , 



s=±i 

-anr>(m, t) + on^im) 



(14) 



In this case, the iteration of mean field equations yield 
results (see Table II) that are much closer to those found 
from the simulations. 

The evolution equations, which we have written as dif- 
ference equations, arc of course nonlinear. In the simplest 
case of asexual reproduction (Eqs. (|[ ) this second or- 
der nonlinearity comes purely from the condition of a 
fixed finite population, and appears in the source term 
for the restoration of the population to its fixed value by 
randomly sampling the asexual population and cloning 
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TABLE II: The distribution of the population with respect 
to number of expressed mutations, obtained from an iteration 
of the mean field equations for Model B. 
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12.7 


Total 
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49.3 



it. With the introduction of sex, the source term in the 
equations for the asexual organisms pi] , |l3| ) acquires a 
contribution from the number of sexual individuals that 
are removed either through death or through sexual re- 
production. Such terms contain nonlinearities up to third 
order. We expect to find nontrivial behaviour in the limit 
of large nonlinearities in these equations, and this turns 
out to be the case, as we explore in the next section. 

IV. LIMITS OF STRONG AND EXTREMELY 
WEAK DRIVING, CHAOTIC BEHAVIOUR 

A. The limit of strong driving 

After the discussion of the last section, it is natural 
to expect that the nonlinearities present in the problem 
should drive it to chaotic behaviour when their amplitude 
is sufficiently large. 

We have tested the limit of V — 1 and found that for 
Model A with hereditary sex, the system becomes unsta- 
ble. The total asexual population and sexual population 
display oscillations with a period of 2 time steps. The 
m-distributions also oscillate for both populations, with 
the same period, the amplitude of the oscillations being 
much larger for the asexuals. For such large values of T, 
at each time step a large number of asexuals are driven to 
large m values and are converted to sexuals, they mate, 
and reduce their expressed mutations. This leads to a 
macroscopic fluctuation in the number of sexuals, with 
the halving of the mating population, which then causes 
a very large number of asexuals to be cloned in turn. 
The time average of the sexual population is depressed 
slightly below the saturation value as a result, as can be 
seen in Fig. (3) . These oscillations are not observed in the 
iteration of the mean field equations. 

A much more striking behaviour is found in Model B 
for large values of a. As we increase the value of a, the 
probability of random conversion to sex, beyond about 
0.05, a spectacular transition takes place to a strange at- 
tractor for the dynamics of both the asexual and sexual 
populations. In place of the well converged m distribu- 
tions for both asexual and sexual populations, shown in 
Figs. 6 one observes that both distributions are inter- 
mittently switching between several meta-distributions. 



4.0 
3.5 




0.5 



°'5.1x10 5 5.15x10' 5.2x1 Q 

t 

FIG. 8: In Model B, we find intermittent variation with time 
of a) < m > a , the average number of mutations for the asex- 
ual population, and b) < m > s , the average number of ex- 
pressed mutations for the sexual population, for the relatively 
large value of a = 0.5. The averages are taken over the pop- 
ulation at time t. F — 10 -3 . It is clearly seen that there are 
two metastable states. The pictures show a window of 10 4 
time steps after the transients are dropped. 



The average value of m computed over the asexual and 
the sexual populations is shown in Fig. 8, and displays 
this striking intermittent behaviour, where the distribu- 
tion of the two populations becomes much more closely 
coupled than in the lower a values. They now move more 
or less in phase, and their excursions take them all the 
way down to the wild type. Now it is only possible to 
talk about a distribution of distributions. To display this 
graphically, we have plotted the distribution of the av- 
erage number of expressed mutations in the two popula- 
tions, (m) a and (m) s , as a function of a. In Fig. 9, we 
show three dimensional plots for these distributions, com- 
piled over 10 4 time steps for each value of a. In Fig. 10, 
a contour plot of the same distribution as in Fig. 9 are 
shown. It is possible to read off from the contour plots 
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that the transition is taking place around a c ~ 0.05. 





FIG. 9: A 3D plot showing the branching distributions of a) 
< m > a b) < m > s with respect to a. After a threshold at 
a ~ 0.05, the distribution displays more than one peak. The 
z-axis indicates the relative weights of these peaks. The total 
population is 1000 and the figure represents a single run of 
10 4 steps. 
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FIG. 10: Contour plot showing the branching of a) < m > a b) 
< m > 3 , as a increases, for a population of 1000, computed 
over 10 4 time steps. 



branch which in the simulations has the smaller weight, 
while the asexuals evolve along the higher (large m) 
branch, which has the greater weight in the simulations, 
and the evolution is completely stable. For a = 0.9 and 
T = 0.1, (m) a = 2.43 and (m) s = 0.47. 



Besides being intermittent, this transition has a dra- 
matic effect on the to distribution of the sexual popula- 
tion, in that it shifts it to much higher values. It can 
be seen in Fig. 10(b) that for a < <r c , the mean m for 
the sexual population is (m) s ~ 0.75, while for large a it 
is comparable to the corresponding value for the asexual 
population, closer to 3. The reason seems to be that with 
the great depletion of the population when too many in- 
dividuals are being switched on to sex and engaging in 
sexual reproduction, the asexuals are cloning too many 
identical copies to make up for the deficit. When these 
are subsequently turned sexual and mate among each 
other, "inbreeding" takes place - there is not sufficient 
genetic diversity for sex to lead to sufficient mixing and 
therefore an amelioration of the effective fitness. 

We have iterated the mean field equations (13|,|l4|) for 
Model B and found that this intermittent behaviour is 
suppressed. The sexuals simply evolve along the lower 



B. The limit of infinitely slow driving (r 

a -> 0) 



or 



In the limit of infinitely slow driving, i.e., r — ► or 
(7 — > 0, we observe a transition to a different phase. 

For r < 1/N, we find a qualitatively different asex- 
ual steady state, where the m distribution has shifted to 
lower m values (compare with Fig. 2 of [Q) and no longer 
has the characteristic minimally stable sand-pile like jlCfl 
distribution. For T — 10~ 4 = (10iV) _1 , over a run of 
10 6 steps, we find n H (m)/N ~ 0.03,0.14,0.44,0.39 for 
to = 0, . . . , 3 respectively, where the peak has moved to 
to = 2 from to = 3, or broadened towards the left. This 
does not seem simply to be due to a slowing down of the 
dynamics. Rather, once the mutation rate drops below 
1/N, the flow over the m = 4 threshold which stabilizes 
the skewed distribution slows down to a dribble. This 
gives the m distribution time to get stabilized at to = 2 
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successively to to = 1,2 and finally to to = 3 where it sta- 
bilizes. Comparison with the simulation results seem to 
indicate that the simulations get stuck at an intermedi- 
ate "metastable" state, while the peak is around m = 2. 
The fact that in the simulation one has to wait around 
until, with a very low probability, a discrete individual is 
pushed over the to = 4 barrier, dies, and is cloned from 
among the live bacteria, while in the mean-field equa- 
tions, there is a weak but steady seepage, which prevents 
this phase transition from taking place. 




10 15 
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FIG. 12: The iterated solutions of the equations for the purely 
asexual population, without the introduction of sex, as a func- 
tion of time for different values of m. V = 10 -4 . 



FIG. 11: The distribution over m for a) asexual b) sexual 
populations, for different values of F for Model A, without 
hereditary sex. The steady state distribution changes and 
the peak on the distribution shifts to a smaller m value as 
one lowers the T value below the threshold 1/N = 10~ 3 . 



rather than being pushed to the m = 3 limit. The mecha- 
nism for the stabilization is provided by the dead bacteria 
being replenished from among the most prevelant extant 
ones. 

Once sex is turned on in Model A, we similarly observe 
that the peaks in the distribution of the asexual and sex- 
ual populations have shifted to lower m values (m = 1 
and to = 2 respectively), as shown in Fig.ll. Although 
the total sexual population is relatively small here, we 
have checked that the fluctuations in the histogram over 
10 different realizations stay small. 

Iteration of the dynamical equations, on the other 
hand, reveal no such phase transition and, for the asexual 
steady state, converge to the same steady state distribu- 
tions as found for T > 1/N. In Fig. 12, we show the time 
series for njj(m) (N = 100) for the asexual population 
without conversion to sex. At time t = 0, the largest den- 
sity is of course at to = 0, and then the maximum shifts 



V. DISCUSSION 

The mechanism of random conversion to sex, in the 
presence of a constant rate of mutations, as investigated 
in this paper as scenario for the maintanence of a macro- 
scopic sexual production, is in fact very closely related 
to "coevolution of cell senescence and diploid sexual re- 
production in unicellular organisms," studied by Cui et 
al. [jy]. In this paper a "senescence clock" ticks off a finite 
lifetime for each bit-string. Sexual reproduction (conju- 
gation) resets the senescence clock; unless this happens 
after a number of generations of cloning, the offspring 
stop dividing and die. 

Our Model B can be seen as a simpler version of the 
model proposed by Cui et al., with an intrinsic mech- 
anism, provided by Muller's ratchet |L2]], for cell senes- 
cence. The constant mutation rate sets the time scale for 
the survival of any given bitstring, unless it succeeds en- 
gaging in sex, with a given probability (our a) . A survival 
function (Eq.(l)) leads to the elimination of genomes car- 
ried by haploid individuals multiplying by asexual repro- 
duction, once they have accumulated too many muta- 
tions as a result of prolonged exposure to the constant 
mutation rate [Q, [l^] . 
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Our Model A goes one step further, in that it makes 
the number of mutations (the cell clock) provide the trig- 
gering mechanism for the transition to diploidy and sex. 
It is gratifying to find that this is a more successful strat- 
egy for establishing a sexual population than a constant 
rate of conversion to sex. 

Chopard et al. |l3| have pointed out that care must be 
taken in the investigation of finite populations, amplifiy- 
ing and stabilizing small fluctuations which in the ther- 
modynamic limit would be attenuated to zero. They em- 
phasize the importance of spatial variations which can- 
not be captured by mean field theories. In this paper we 
have demonstrated the relevance for finite populations of 
discrete stochastic events, whose effect in the very weak 
driving limit cannot be captured by the "mean field" 
equations. In the very weak driving limit the system 
is below the hydrodynamic regime, and exhibits a qual- 
itatively different phase than which is described by the 
continuum approximations. 

In a recent article Pekalski [[u] has studied a model 
which is in many ways similar to ours. There the suc- 
cess of sexual reproduction, meiotic parthenogenesis and 
asexual reproduction, in maintaining a finite population 



in the face of periodically changing environmental condi- 
tions and a constant mutation rate, is studied in terms of 
the relative sizes of the populations. Age is included in 
the model as a parameter which reduces the fitness. The 
populations do not interact. The findings are that mei- 
otic parthenogenesis and sexual reproduction are more 
favorable than mitotic reproduction, with slight differ- 
ences between them depending on the precise conditions. 

Further work is in progress, to investigate the effect of 
finite temperature, and of including the possibility of ge- 
netic crossover and meiotic parthenogenesis, in our mod- 
els. Results on the autonomous viability of the sexual 
population, after the steady conversion from the haploid 
population has been switched off (but mitosis allowed for 
the diploids), will be reported in a future publication. 
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